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A  method  for  resolving  inverse  heat  conduction  problems  using  approximating  polynomials  is  described.  The 
coefficients  are  obtained  by  constraining  the  even-numbered  spatial  derivatives  of  the  polynomial  at  two  points 
within  the  domain  using  the  heat  equation.  The  result  is  a  linear  time-invariant  system  that  can  be  implemented  as  a 
digital  filter  and  is  suitable  for  application  in  a  real-time  sensor.  The  method  places  no  restrictions  on  boundary  or 
initial  conditions,  and  temperature  dependence  of  transport  properties  is  accounted  for  in  an  approximate  way.  The 
performance  of  the  model  has  been  validated  using  analytical  and  numerical  solutions,  and  results  are  presented  in 
the  frequency  and  time  domains  for  several  orders  of  polynomials.  A  propagation  of  error  analysis  is  presented  and  is 
used  to  establish  the  optimum  location  of  the  sensors. 


Nomenclature 

Ojj  =  matrix  of  location-dependent  coefficients,  a 13, 
a14:m2;  a2 1,  a22:m -1;  a22,  a^'.vcr,  a43,  a44:m_1 
bj(t)  =  vector  of  temperature  measurements,  b{,  b2:K;  b3, 
b4:  K/m2 

eft)  =  vector  of  temperature  profile  polynomial  coefficients, 
cqiK;  c2:  K/m;  c3:K/m2;  c4:K/m3 
cp  =  heat  capacity,  J/kg  •  K 

Im  =  imaginary  part 

i  =  imaginary  number,  x/— T 

k  -  thermal  conductivity,  W/m  •  K 

Pn  =  temperature  profile  polynomial  of  degree  n,  K 

q  =  heat  flux,  W/m2 

Re  =  real  part 

T  =  temperature,  K 

T  =  rate  of  change  of  temperature,  K/s 
t  =  time,  s 

x  =  distance  from  the  surface,  m 

a  =  thermal  diffusivity,  m2/s 

=  low-pass  filter  coefficient  for  time  derivative,  s~* 

T  =  gain,  nondimensional 

Yj  =  low-pass  filter  coefficient  for  time  average 

At  =  time  step,  s 

S  =  uncertainty 

p  =  density,  kg/m3 

a  =  length  scale  parameter,  nu1 

(p  =  phase  angle,  rad 

co  =  angular  frequency,  rad/ s 

|  =  amplitude 


Subscripts 

i,  j  =  polynomial  and  matrix  indices 

m  =  low-pass  filter  index 

n  =  order  of  polynomial  approximation  for  temperature 

1  =  position  nearest  to  the  measurement  surface 

2  =  position  farthest  from  the  measurement  surface 
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%  =  nondimensional  number 


Superscripts 

=  frequency  domain  quantity 
=  low-pass  filtered  quantity 
*  =  complex  conjugate 


I.  Introduction 

THE  measurement  of  the  transient  heat  flux  and  the  surface 
temperature  in  heat-sink  combustion  chambers  continues  to 
present  technical  challenges  to  the  instrumentation  engineer.  Sensor 
failure  rates  are  high,  and  measurement  accuracies  and  uncertainties 
are  not  well  characterized.  These  shortcomings  have  had  a  significant 
impact  on  some  recent  programs  that  have  used  heat-sink  test  articles 
to  acquire  data  for  the  validation  of  heat  transfer  predictions  at  liquid 
rocket  engine  operating  conditions  [J_,2]. 

There  are  numerous  types  of  heat  flux  sensors,  but  a  relatively 
small  subset  is  capable  of  operating  in  rocket  chamber  conditions  in 
which  the  heat  flux  levels  can  exceed  10s  W/m2,  and  the  surface 
temperatures  of  1000  K  are  typical.  Diller  [3]  reviewed  the  devices 
that  have  been  used  and  organized  them  into  methods  that  relied  on 
temperature  differences  over  a  spatial  distance  with  known  thermal 
resistance  and  temperature  differences  over  time  having  known 
thermal  capacitance.  The  most  commonly  used  method  has  been  the 
coaxial  thermocouple,  which  is  an  example  of  the  second  type.  A 
thermocouple  junction  is  formed  on  the  surface  of  the  chamber 
between  a  wire  of  one  type  of  thermocouple  material  and  a  sur¬ 
rounding  sheath  of  another  type.  The  heat  flux  is  determined  from  the 
measured  temperature  boundary  condition  using  a  one-dimensional 
transient  solution  to  the  heat  equation.  The  junction  is  typically  a  few 
micrometers  in  thickness  and  is  often  formed  by  lightly  scratching 
the  surface  to  drag  filaments  of  one  type  of  material  across  the  elec¬ 
trically  insulating  layer  to  the  other  type.  In  some  applications,  when 
erosion  of  the  surface  occurs,  the  junction  is  continuously  reformed, 
and  this  has  led  to  the  description  of  coaxial  thermocouples  as 
eroding  thermocouples.  However,  in  heat-sink  chambers,  it  is  quite 
common  to  find  that  the  junction  disappears  at  some  point  during  a 
test,  and  the  sensor  fails. 

Other  methods  have  been  developed  that  do  not  rely  on  surface 
temperature  measurements  but  embed  the  sensors  within  the  wall 
where  they  are  protected  from  erosion.  In  the  null-point  calorimeter,  a 
hole  is  drilled  from  the  backside  of  the  chamber  wall  and  a  ther¬ 
mocouple  is  inserted.  The  bead  is  brazed  or  resistance  welded  to  the 
bottom  of  the  hole.  Null  point  refers  to  a  distance  from  the  bottom  of 
the  hole  to  the  inner  wall,  where  the  disturbance  to  the  flow  of  heat 
caused  by  the  hole  results  in  the  junction  reading  nearly  equal  to  the 
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inner  wall  temperature.  The  construction  of  null-point  calorimeters  is 
challenging.  The  junction  cannot  be  visually  inspected,  and  large 
measurement  errors  can  result  from  manufacturing  flaws  [4], 

Another  method  using  embedded  temperature  sensors  is  the  plug- 
type  heat  flux  gauge  of  Liebert  [5],  An  annular  groove  is  machined 
into  the  chamber  wall  to  form  a  post,  and  thermocouples  are  attached 
at  several  axial  points  along  the  outside  of  the  post.  A  polynomial 
curve  is  used  to  extrapolate  the  temperatures  to  the  wall  position,  and 
an  integral  method  is  used  to  calculate  the  total  heat  load  to  the  plug 
from  transient  temperature  measurements.  Two-dimensional  effects 
can  be  significant  in  this  type  of  device.  The  dimensions  of  the  groove 
are  critical,  and  significant  errors  can  result  from  the  disturbance  to 
the  flow  of  heat  [6]. 

Recently,  Conley  et  al.  used  embedded  temperature  sensors  to 
measure  transient  heat  flux  in  a  subscale  combustion  chamber  [7], 
Transient  temperature  data  were  reduced  to  heat  flux  using  an 
equation  developed  from  an  energy  balance  on  a  control  volume 
between  the  sensors.  The  equation  includes  terms  for  Fourier’s  law  of 
heat  conduction  and  heat  storage.  Nominal  constant  values  were  used 
for  the  properties. 

The  data  reduction  methods  developed  by  Leibert  [5]  and  Conley 
et  al.  [7]  are  particular  examples  of  a  more  general  class  of  inverse 
solution  procedures.  Inverse  heat  transfer  problems  are  challenging 
because  surface  conditions  must  be  obtained  from  temperature  sen¬ 
sors  embedded  within  objects  that  experience  attenuated  and  time- 
lagged  responses  to  changes  in  the  boundary  conditions.  Solution 
methods  must  be  stable  and  accurate  in  the  presence  of  noisy  signals. 
Most  inverse  methods  are  based  on  minimizing  the  sum  of  square 
residuals  (SSR)  of  the  difference  between  measurements  and  a 
solution  of  the  heat  equation.  An  initial  estimate  is  made  for  the 
boundary  condition,  and  the  response  at  the  sensor  location  is  cal¬ 
culated.  Then  the  estimate  is  improved  using  a  gradient  method  until 
a  convergence  criterion  for  SSR  is  satisfied.  The  method  used  to 
calculate  the  response  is  specific  to  the  problem  and  can  be  any  of  the 
analytical  and  numerical  methods  that  exist  for  solving  the  heat 
equation.  However,  there  is  a  fundamental  physical  limit  on  the  size 
of  the  time  step  because  there  must  be  sufficient  time  for  the  change  in 
the  boundary  to  propagate  to  the  sensor  location.  A  primary  focus  of 
research  in  this  area  is  the  development  of  methods  that  can  achieve 
the  smallest  possible  time  step  without  incurring  instability.  A  widely 
used  procedure  is  the  function  specification  method  [8].  A  functional 
form  for  the  time  dependence  of  the  boundary  condition  is  assumed, 
and  the  parameters  of  the  function  are  then  determined  by  mini¬ 
mizing  SSR.  The  function  may  span  the  entire  time  domain  or  it  may 
be  parsed  into  subintervals.  Another  widely  used  class  of  technique  is 
based  on  regularization  [8],  These  techniques  add  additional  factors 
to  the  least  squares  function  to  enhance  stability  and  convergence. 

If  the  method  used  to  solve  the  heat  equation  is  a  numerical 
technique,  the  number  of  processor  operations  will  likely  require  that 
the  inverse  calculations  be  performed  as  a  postprocessing  operation. 
However,  if  the  method  is  based  on  an  analytical  solution,  it  may  be 
possible  to  perform  the  calculation  in  real  time.  Inverse  methods 
based  on  analytical  solutions  are  generally  limited  to  situations  for 
which  thermal  properties  are  effectively  constant  and  specific  bound¬ 
ary  conditions  are  satisfied.  An  important  subset  of  these  methods 
has  been  developed  for  the  case  in  which  the  transient  surface 
temperature  is  measured  and  the  test  duration  is  short,  such  that  the 
body  can  be  treated  as  an  infinitely  thick  slab.  The  method  of  Cook 
and  Felderman  assumes  that  heat  flux  is  constant  within  a  sampled 
time  interval  and  the  analytical  relation  between  heat  flux  and  surface 
temperature  is  used  to  obtain  the  surface  flux  from  the  measured 
temperatures  as  a  summation  of  contributions  from  each  time  interval 
[9].  Cook  also  suggested  a  modification  to  the  method  to  account  for 
temperature-dependent  properties  [10], 

There  have  been  a  number  of  publications  describing  the  applic¬ 
ation  of  digital  signal  processing  techniques  to  the  measurement  of 
heat  transfer.  Beck  et  al.  discusses  the  use  of  an  efficient  digital  filter 
algorithm  for  linear  problems  and  its  application  to  online  or  real¬ 
time  data  processing  [8],  In  the  method  by  Beck  et  al.,  the  filter 
coefficients  can  be  derived  for  the  particular  problem  from  any  of 
the  numerical  or  analytical  approaches.  Once  obtained,  the  filter 


coefficients  do  not  change,  and  the  heat  flux  calculation  is  reduced  to 
the  product  of  two  arrays,  although  a  specific  origin  of  time  is 
required.  Marineau  and  Homung  described  an  efficient  method  for 
calculating  heat  flux  from  surface  temperature  measurements  based 
on  digital  signal  processing  techniques  [11],  The  relationship  be¬ 
tween  heat  flux  and  surface  temperature  is  based  on  the  Green’s 
function  for  the  impulse  problem,  and  the  heat  flux  is  deconvolved 
from  the  Green’s  function  and  the  surface  temperature  in  the 
frequency  domain.  Oldfield  developed  an  efficient  method  for 
processing  signals  from  thin  film  heat  flux  gauges  based  on  the  linear 
time-invariant  (LTI)  system  theory  [121.  The  analytical  solution  for 
the  impulse  response  of  the  sensor  is  convolved  with  the  measured 
temperature  distribution  to  obtain  the  surface  heat  flux. 

Frankel  et  al.  have  suggested  other  approaches  that  could  be  used 
for  real-time  processing  [13].  For  a  one-dimensional  problem,  a 
Taylor  series  can  be  used  to  extrapolate  from  the  measurement  point 
to  the  surface  condition,  with  the  spatial  derivatives  in  the  series 
obtained  from  rate  sensors  for  temperature  and  heat  flux  (which  is 
nonintrusive  to  the  flow  of  heat)  and  the  heat  equation.  Frankel  also 
suggested  a  two-sensor  arrangement ,  wherein  a  Taylor  series  for  the 
first  sensor  is  expanded  about  location  x,  and  a  series  for  the  second 
sensor  is  expanded  about  2x.  Adding  the  two  expansions  eliminates 
terms  and  results  in  an  order  of  accuracy  equal  to  the  one  sensor 
expansion,  with  derivatives  in  time  one  order  lower. 

In  this  paper,  we  extend  these  ideas  and  report  a  number  of  new 
results  for  the  class  of  methods  based  on  the  extrapolation  of  a 
function  for  the  two-sensor  arrangement.  We  have  found  that  this 
method  leads  to  very  satisfactory  solutions  to  the  major  technical 
challenges  of  accuracy,  time  response,  stability  in  the  presence 
of  noise,  uncertainty,  computational  efficiency,  failure  rate,  and 
producibility  in  the  heat-sink  chambers  used  for  rocket  heat  transfer 
experiments. 


II.  Derivation  of  Model 

The  domain  of  interest  is  a  continuous  one-dimensional  body  with 
smoothly  varying  or  constant  thermal  properties  (Fig.  j_).  The  initial 
temperature  distribution  is  assumed  to  be  continuous  and  smoothly 
varying.  The  boundary  conditions  at  the  surfaces  are  arbitrary  and 
may  be  discontinuous  in  time.  The  principle  goal  of  the  analysis  is  to 
determine  the  transient  temperature  and  heat  flux  at  x  =  0  from 
temperature  data  at  xt  and  x2. 

We  begin  by  assuming  a  one-dimensional  flow  of  heat  in  the  x 
direction  and  approximate  the  temperature  profile  in  the  body  with  a 
polynomial  in  x  with  time-dependent  coefficients: 

n  + 1 

Pn(x,  t)  =  'Y^ci(t)xi~l  (1) 

;= l 

In  addition  to  matching  temperatures  at  the  two  points,  the  poly¬ 
nomial  is  also  required  to  satisfy  the  one-dimensional  heat  equation. 
Specifically,  it  is  required  to  match  the  second  derivative  of  tem¬ 
perature  with  respect  to  position: 


Fig.  1  Domains  addressed  by  current  theory:  a)  finite  slab,  b)  finite  slab 
with  sensor  on  second  surface,  c)  semi-infinite  slab,  and  d)  internal 
points. 
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d2  T  _  1  dT 

dx 2  a  dt 


(2) 


As  shown  by  Burggraf  [15],  the  heat  equation  can  be  extended  to 
any  even-numbered  derivative  of  x\ 


d2,„T  _  1  dmT 


(3) 


Expressions  for  the  c,-(f)  in  Eq.  (1)  are  obtained  by  equating 
Pn(x ,  t)  and  its  even-numbered  spatial  derivatives  with  measured 
values.  For  the  case  of  P3  (x,  t),  the  set  of  equations  can  be  expressed 
in  matrix  form  as  follows: 


T(xu  t ) 
T(x2,  t) 

s2nxt.,) 

dx2 

S2T(x2,  t) 

L  3x2  -I 


1  X\  x\ 

1  x2  x\ 

0  0  2 

0  0  2 


The  inverse  is  the  following: 


*2  -*1 

xix2(2x2-xi) 

Ci(0 

*2~*1  x2  xl 

6(X2-Xi  ) 

c2(t) 

-1  1 

x*— 2xi  *2  —2*1 

— 

*2-*l  *2“ *1 

6(x2-xi  ) 

C3(t) 

c4(t ) _ 

0  0 

0  0 

x2 

2(*2-*i) 

-1 

6(*2-*i) 

'ci(f) 

c2(l) 

C3(t) 

_  c4(t)  _ 

(4) 


*l*2(*2~2*l) 

6(x2-xi) 

2*i +2*i  *2—  Xj 

6(x2-xi) 

-xi 

2(x2-xi) 

1 

6(x2— *i) 


T(x2,  t) 


Pt(x,j) 

a? 

82r(x2.t) 

L  a?  J 


(5) 


c,(0  =  a,jbj(t) 


(6) 


All  terms  in  djj  depend  only  on  the  sensor  locations  and  need  to  be 
calculated  once.  The  vector  bj(t)  contains  the  experimental  mea¬ 
surements.  The  second  derivative  terms  in  bj(t)  are  obtained  from 
Eq.  (2).  The  method  used  to  calculate  the  time  derivative  in  Eq.  (2)  is 
critical  to  the  success  of  the  technique  because  differentiation  of 
experimental  data  can  be  a  highly  unstable  process  [13, 14, _16,  T7]  and 
low-pass  filtering  is  usually  required.  The  vector  c,(t)  can  be  used  to 
reconstruct  an  approximation  for  the  temperature  profile  in  the 
sensor,  but  the  most  important  application  is  in  the  estimate  of  the 
boundary  conditions: 

T(0,t)*tP3(f),t)  =  cl(t)  (7) 

9(0,  t)  ^  -fc8^0’  f>  =  -kc2(t)  (8) 

dx 


Fig.  2  Polynomial  approximations  to  temperature  profile  in  a  copper 
slab  exposed  to  steady-periodic  heat  flux  at  left  boundary.  Sensors  are 
located  at  2  and  5  mm. 


second  and  fourth  derivatives  and  is  indistinguishable  from  the  exact 
solution.  The  maximum  fractional  error  at  the  surface  is  0.03. 

A.  Frequency  Domain  Analysis 

The  equations  in  (5)  constitute  an  LTI  system.  The  polynomial 
approximation  is  obtained  from  linear  combinations  of  the  mea¬ 
surements  contained  in  the  bj(  t)  vector  with  the  constant  coefficients 
of  the  djj  matrix.  The  system  is  time  invariant  because  the  output  is 
based  on  current  values  of  the  measurements  only.  Therefore,  a  time 
shift  in  the  input  produces  an  equal  time  shift  in  the  output.  LTI 
systems  can  be  completely  represented  by  an  impulse  response  or 
equivalently,  in  the  frequency  domain,  by  the  gain  and  phase  re¬ 
sponse.  In  this  section,  we  present  the  frequency  domain  behavior  for 
approximating  polynomials  of  orders  1,  3,  and  5.  To  avoid  lengthy 
expressions,  we  illustrate  the  transformation  using  just  the  P3  ap¬ 
proximation.  We  begin  with  the  following  expression  for  the  time- 
domain  surface  temperature,  which  has  been  obtained  by  replacing 
the  second  derivatives  with  respect  to  x  in  the  bj(t)  vector  with  the 
right-hand  side  of  Eq.  (2): 

T(xut)  t(x2,  t) 

P3( 0,  t)  =  anT(xu  t)  +  dvT(x2,  t)  +  al3 - h  al4 - - — 

a  a 

(9) 

Following  Cole,  we  transform  to  the  frequency  domain  using  the 
complex-valued  T(x,co)  [181: 

T(x,  t)  =  Re[f(x ,  co)ei,M]  (10) 


Note  that  there  are  no  specific  limitations  that  have  been  placed  on 
the  boundary  conditions.  The  model  can  be  applied  to  semi-infinite 
or  finite  thickness  slabs.  The  sensor  at  x2  can  be  embedded  in  the 
body  or  it  can  be  placed  on  the  surface.  The  surface  can  be  adiabatic 
or  actively  cooled  at  a  known  or  unknown  rate.  The  technique  can  be 
used  to  estimate  temperature  and  heat  flux  profiles  within  a  body. 
The  extension  of  this  approach  to  higher-order  polynomials,  using 
Eq.  (3),  and  radial  geometry  is  straightforward. 

As  an  initial  illustration  of  the  ability  of  the  polynomial  approxi¬ 
mations  to  reproduce  temperature  profiles  and  boundary  conditions, 
we  take  as  an  example  a  uniform  slab  of  copper  exposed  to  a  steady- 
periodic  heat  flux  boundary  condition.  The  heat  flux  is  oscillating 
sinusoidally  at  10  Hz,  and  temperature  sensors  are  located  at  depths 
of  2  and  5  mm  from  the  heated  surface.  For  the  purpose  of  this 
illustration,  the  sensors  are  free  of  noise.  In  Fig.  2,  the  analytical 
solution  is  compared  with  polynomial  approximations  of  orders  1 ,  3 
and  5.  The  order  1  approximation  is  a  simple  lineal'  extrapolation  of 
temperature  to  the  surface.  It  is,  in  effect,  a  steady  approximation  and 
exhibits  a  large  error  and  time  delay.  The  maximum  fractional  error  at 
the  surface,  based  on  the  amplitude  of  the  input  wave,  is  0.81.  The 
order  3  approximation  includes  the  second  derivative  terms,  and  the 
results  are  much  improved.  The  maximum  fractional  error  at  the 
surface  is  reduced  to  0.11.  The  order  5  approximation  includes 


Substituting  Eq.  (10)  into  Eq.  (9), 


P3(0,  co)  =  ^a„  +  t^ai3jf(xi,  co)  +  ^a12  +  i°\l4^T(x2,  co) 


(ID 


The  corresponding  expression  for  surface  heat  flux  is 


—  k 


dP3( 0,  co) 
dx 


=  —  k  (  a2 ,  +  i  —  a23  )7Yxi,  «) 

LV  "  «  J 


(a22  +  i°^a24^T(x2,  co) 


(12) 


To  illustrate  the  behavior  of  the  model,  we  obtain  inputs  from  an 
exact  analytical  solution.  We  consider  the  case  of  a  uniform  semi¬ 
infinite  slab  with  a  steady-periodic  heat  flux  boundary  condition.  The 
solution  to  this  problem  is  [19]: 

\a\e~ax 

T(x,  co)  = - ,  (0  <  x  <  oo,  co  >  0)  (13) 

k  a 


q(x,  co)  =  \q\e 


(0  £  x  <  oo,  co  >  0) 


(14) 
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tr=  (1  +  0 


(15) 


Gain  is  defined  as  the  ratio  of  the  amplitude  of  the  model  to  the  actual 
value: 


r  (Pn)  =  \P„(0,  «)|/|T(0,  w)|  (16) 

Phase  is  defined  as  the  difference  between  the  phase  of  the  model  and 
the  actual  value: 


<t>{Pn)  =  tan  1 


Im[-P„(0,  &>)]\ 
Re[Pn(0,  a>)]) 


-tan-'(Imlf(°’W)]) 
\Re[T( 0,  co)]) 


(17) 


The  gain  and  phase  expressions  for  surface  heat  flux  are  identical  to 
those  for  surface  temperature. 

Gain  and  phase  plots  are  shown  in  Figs.  3-6  for  polynomials  of 
orders  1,  3,  and  5.  Frequency  has  been  nondimensionalized  using  a 
time  scale  based  on  the  thermal  diffusivity  and  the  depth  of  the  near¬ 
surface  sensor: 


The  value  of  x2/x\  in  each  case  is  2.3.  This  value  was  determined 
to  be  near  optimal  and  is  discussed  further,  next,  in  the  context  of  the 
noise  sensitivity.  The  gain  behavior  for  surface  temperature  is  shown 
in  Fig.  3.  All  orders  of  approximation  converge  to  a  gain  of  one  at  low 
frequency  but  deviate  as  frequency  is  increased.  The  Pt  approxi¬ 
mation  rolls  off  slowly  to  zero  and  has  a  —3  db  frequency  of  2.  The 
P3  approximation  exhibits  a  region  of  gain  greater  than  1  and  then 
rolls  off  at  a  frequency  of  20.  The  P5  approximation  has  a  larger  gain 
region  and  a  — 3  dB  frequency  of  1 00.  The  region  of  gain  greater  than 
one  can  be  eliminated  by  low-pass  filtering,  as  will  be  shown  later. 
The  gain  at  high  frequency  would  manifest  itself  in  the  time  domain 
as  an  overshoot  in  the  response  when  subjected  to  rapid  changes  in 


Fig.  4  Surface  temperature  phase  for  three  orders  of  polynomial 
models. 


Fig.  5  Surface  heat  flux  gain  for  three  orders  of  polynomial  models. 


Fig.  6  Surface  heat  flux  phase  for  three  orders  of  polynomial  models. 


the  inputs.  The  phase  behaviors  show  similar  points  of  departure,  but 
all  approximations  exhibit  monotonically  increasing  phase  lag 
(Fig.  4).  From  these  plots,  we  can  state  that  the  P3  approximation  has 
a  frequency  response  approximately  one  order  of  magnitude  higher 
than  the  P ,  approximation,  and  the  P5  approximation  is  approxi¬ 
mately  two  orders  higher. 

The  behavior  of  the  surface  heat  flux  plot  (Fig.  5)  is  qualitatively 
similar  to  surface  temperature,  but  the  points  of  departure  and  roll  off 
are  shifted  to  lower  frequencies.  The  —  3  db  point  of  the  P{  ap¬ 
proximation  occurs  at  a  frequency  of  0.3  because  the  simple  linear 
extrapolation  does  a  poor  job  capturing  the  slope  of  the  temperature 
profile  at  the  surface.  The  P3  approximation  shows  smaller  amplitude 
in  the  gain  region  than  the  surface  temperature  plot,  but  the  P5 
approximation  continues  to  have  a  region  of  large  gain.  The  phase 
characteristics  of  heat  flux  are  also  monotonically  decreasing 
(Fig.  6).  The  point  of  departure  for  P{  is  shifted  to  a  frequency  of  less 
than  0.01 .  The  P3  approximation  begins  to  exhibit  a  phase  lag  at  a 
frequency  of  0.3,  but  the  P5  does  not  exhibit  a  phase  lag  until  20. 
From  these  plots,  we  can  state  that  the  P3  approximation  has  a 
frequency  response  at  least  three  orders  of  magnitude  higher  than  the 
P{  approximation,  and  the  P5  approximation  is  approximately 
four  orders  higher.  These  improvements  come  about  solely  as  a  result 
of  a  higher-order  polynomial  approximation  to  the  temperature 
profile. 


B.  Low-Pass  Filtering 

The  formulation  as  an  LTI  system  facilitates  the  use  of  digital 
signal  processing  techniques.  Low-pass  filtering  of  the  input  signals 
can  improve  the  performance  of  the  polynomial  models  by  sup¬ 
pressing  the  propagation  of  errors  from  the  input  signals  and  by 
eliminating  the  region  of  gain  greater  than  one.  The  design  of 
optimum  filters  is  an  extensive  subject  and  has  been  treated  by  several 
authors  in  the  context  of  inverse  heat  conduction  methods.  Hensel 
showed  that  filtering  the  input  data  improves  the  stability  of  an 
inverse  calculation,  and  he  recommended  the  use  of  a  Gaussian  filter 
because  it  does  not  distort  phase  and  has  an  acceptable  rate  of  roll  off 
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Fig.  7  Heat  flux  gain  for  P5  model  with  low-pass  filter  for  four  values  of 
filter  window  width. 


Fig.  8  Heat  flux  phase  for  P5  model  with  low-pass  filter  for  four  values 
of  filter  window  width. 


at  high  frequency  [20].  Frankel  used  a  Gaussian  filter  to  suppress 
noise  in  the  input  signals  before  using  a  finite  difference  calculation 
to  obtain  a  time  derivative  [141.  Frankel  suggested  a  method  for 
interrogating  the  power  spectrum  to  locate  an  optimum  cutoff  fre¬ 
quency  for  the  input  filter.  However,  there  remained  a  significant 
propagation  of  error  to  the  derivative,  which  he  attributed  to  the  ill- 
posed  nature  of  numerically  differentiating  noisy  data  [14], 

An  important  factor  in  choosing  a  low-pass  filter,  in  addition  to  the 
phase  and  roll-off  characteristics,  is  the  efficiency  with  which  the 
calculation  can  be  performed.  The  most  efficient  methods  are  im¬ 
plemented  as  convolutions  in  the  time  domain: 

m 

T(xn  tn)  =  Y  YjT(Xi,  tn+J)  (19) 

J=-m 


to  obtain  the  following  forms  for  the  surface  temperature  and  heat 
flux: 


^(O,  tK)  = 


tn+J ) 


m  r  /  \ 

£  [a"yj  +  C'a  NT(x" 

j=-m L  V  / 

+  ^ai2 Yj  +  ~ Pj^T(x2,  tn+j) 

dP3( 0,  tn)  \(  a23  a  \rr,  ,  , 

k - Yx -  =  ~kl^  lfl21  Yj+—Pj)T(Xl’  tn+j) 

j=-m L  V  / 

+  (a22 Yj  +  tn+j) 


(21) 


(22) 


_  m 

t(x„  tn)  =  Y  PiT(Xu  tn+j )  (20) 

j=-m 

A  filter  that  is  well  suited  for  this  type  of  processing  is  the 
polynomial  smoothing  filter,  also  known  as  the  Savitzky-Golay  filter 
[2J..22] .  The  linear  combination  of  the  coefficients  Yj  and  /3  j ,  with  the 
data  T(xh  f„+J),  results  in  filtered  values  that  are  identical  to  those 
obtained  from  a  least  squares  fit  of  a  polynomial.  However,  the  filter 
is  computationally  efficient  because  it  eliminates  the  need  for  a 
matrix  inversion. 

It  is  informative  to  transform  the  low-pass  filter  into  the  frequency 
domain  in  order  to  characterize  the  effect  on  the  gain  and  phase 
behavior  of  the  polynomial  model.  This  is  accomplished  by  applying 
the  discrete  Fourier  transfonn  to  Eqs.  ( _19)  and  (20)  and  then  using  the 
filtered  frequency  domain  temperatures  in  Eqs.  (11)  and  (12)  [23]. 
The  effect  of  a  quadratic  polynomial  filter  on  the  heat  flux  response  is 
shown  in  Figs.  7  and  8  for  the  P5  model.  The  figures  contain  curves 
for  four  values  of  the  filter  window  width  based  on  the  characteristic 
time  scale  for  the  diffusion  of  heat  x\/a.  The  window  width  of 
0.25a/ x\  is  too  narrow  to  remove  the  high  frequencies  in  which  the 
gain  is  greater  than  one.  The  window  width  of  0.52a /x\  removes  the 
region  of  gain  greater  than  one  and  approaches  most  closely  to  the 
ideal  gain  curve.  Further  increasing  the  window  width  to  x\/a 
reduces  the  frequency  response  by  approximately  a  factor  of  3.  The 
filter  introduces  characteristic  high  frequency  modes,  and  these  are 
sometimes  cited  as  a  drawback  of  the  Savitzky-Golay  filter.  How¬ 
ever,  for  this  application,  although  the  modes  create  large  errors  in 
phase  above  a  frequency  of  20,  the  gain  is  negligible  and  the  results 
are  not  significantly  affected. 

C.  Efficiency 

The  main  advantage  of  the  current  method  over  the  more  general 
inverse  methods  is  that  it  can  be  formulated  as  a  digital  filter  and 
calculated  in  the  time  domain  with  relatively  few  processor  opera¬ 
tions.  To  obtain  the  minimum  number  of  operations,  the  filter 
coefficients  can  be  combined  with  the  atj  coefficients.  Using  the  P3 
case  as  an  example,  Eqs.  (19)  and  (20)  can  be  substituted  into  Eq.  (5) 


The  coefficients  multiplying  the  temperatures  can  be  calculated 
once,  and  then  only  4 m  +  2  multiplications  and  2m  +  1  additions 
are  required  to  evaluate  surface  temperature  and  heat  flux  (or  about 
100  processor  operations  for  a  typical  value  of  m  ss  10).  As  an 
illustration.  Figs.  9-1 1  contain  coefficients  for  31  point  filters  for 
surface  heat  flux  for  polynomials  of  orders  1,  3  and  5,  respectively. 
Note  also  that  the  number  of  operations  is  not  affected  by  the  order  of 
the  approximating  polynomial.  The  order  affects  only  the  value  of  the 
coefficients. 


D.  Noise  Sensitivity 

The  propagation  of  noise  from  the  temperature  measurements  to 
the  predicted  surface  temperature  and  heat  flux  can  be  evaluated  by 
treating  each  temperature  measurement  as  an  independent  random 
variable  and  using  the  square-root-of-sum-of-squares  approach  for 
error  propagation  [24].  In  the  following,  we  assume  that  the 
uncertainty  in  temperature  measurement  ST  is  the  same  at  the  two 
measurement  locations,  and  the  uncertainties  in  measurement  posi¬ 
tions  and  material  properties  are  negligible.  Proceeding  from 
Eqs.  (21)  and  (22)  for  the  P}  approximation,  we  have 


^3(0,  tj  = 

+ 


m  /  \  2 

E(a»Jo+T^) 872 

\=—m  y  / 

(°12  Yj+^-Pj)872 
WJ)  =  KP3(Pt  Q][yJo*Yj  + 

+  («22 Yj+~Pj^  &1 


(23) 


(24) 


Equations  (23)  and  (24),  as  well  as  the  equivalent  equations  for  the 
other  orders  of  approximation,  can  be  rewritten  in  terms  of  the 
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j,  time  index 


Fig.  9  Coefficients  for  calculation  of  surface  heat  flux  based  on  Pl 
approximation. 


Fig.  10  Coefficients  for  calculation  of  surface  heat  flux  based  on  P3 
approximation. 


Figures  12  and  _L3  contain  the  results  for  the  three  orders  of 
approximation  plotted  as  functions  of  the  depth  ratio  of  the  sensors 
x2/xl.  The  noise  increases  with  the  order  of  the  model  because  the 
magnitude  of  the  coefficients  increases,  as  can  be  seen  in  Figs.  9-H . 
The  P  |  model  has  the  lowest  level  of  noise  and  it  decreases 
monotonically  with  the  depth  ratio.  The  P3  and  P5  approximations 
both  exhibit  minima  with  respect  to  x2/x\  for  both  surface  tem¬ 
perature  and  heat  flux.  Flowever,  the  minima  are  broad,  and  values  in 
the  range  of  2. 3-2. 8  are  effectively  equivalent.  Noise  increases  for 
x2/ X\  <  2,  and  so  this  range  should  be  avoided.  An  important 
observation  is  that  in  all  cases  the  noise  is  bounded.  The  method  is 
inherently  stable  because  it  has  the  characteristics  of  a  low-pass  filter 
as  discussed  previously. 

E.  Temperature-Dependent  Properties 

The  linear  form  of  the  heat  equation,  given  in  the  preceding  section 
as  Eq.  (2),  is  valid  when  thermal  properties  are  constant.  Flowever,  in 
many  applications  temperature  dependence  is  significant.  The  most 
basic  example  is  the  case  of  a  sensor  that  undergoes  a  large 
temperature  change  during  the  period  of  measurement.  In  heat-sink 
rocket  chamber  experiments,  the  wall  temperature  can  change  by 
800  K  or  more,  resulting  in  a  30%  decrease  in  the  thermal  diffusivity. 
Spatial  gradients  in  thermal  properties  can  also  be  significant.  As  a 
first  approximation,  when  the  temperatures  at  the  two  measurement 
locations  are  significantly  different,  the  local  temperatures  can  be 
used  to  evaluate  the  thermal  diffusivity  and  the  second  derivative 
term  in  Eq.  (2).  An  improved  approximation  including  spatial 
gradients  in  properties  can  be  obtained  through  the  nonlinear  fonn  of 
the  one-dimensional  heat  equation: 


Fig.  11  Coefficients  for  calculation  of  surface  heat  flux  based  on  P5 
approximation. 


Fig.  12  Noise  in  surface  temperature  prediction. 


following  signal-to-noise  ratios,  such  that  the  only  independent 
variables  are  the  depth  ratio  x2/x{  and  the  filter  parameters.  The  filter 
parameters  used  in  producing  the  plots  were  those  for  a  quadratic 
smoothing  function  with  a  window  width  of  Q.52a/x\  and  m  =  10. 
The  particular  values  are  not  significant  because  we  are  primarily 
interested  in  the  relative  levels  of  noise  for  the  three  approximations. 
Noise  levels  on  the  inputs  can  always  be  reduced  by  increasing  the 
sampling  rate  in  order  to  obtain  more  points  within  the  specified  time 
window: 


_  SPn( 0,  t) 
ST 

S(kdpjdx) 


(25) 


kST /x1 


(26) 
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In  this  case,  the  second  derivative  becomes 

g2r  =  i9r_idfc/9ry 

dx2  a  dt  kdT\dx)  ' 

To  a  first  approximation,  the  first  term  on  the  right-hand  side 
increases  proportionally  with  the  heat  flux,  whereas  the  second  term 
increases  with  the  heat  flux  squared.  For  pure  copper,  the  second  term 
becomes  significant  when  heat  flux  is  on  the  order  of  5  x  107  W/m2. 
The  temperature  gradient  in  the  second  term  cannot  be  measured 
directly  but  can  be  estimated  from  the  linear  form  of  the  heat 
equation,  and  this  can  be  used  in  Eq.  (28)  to  obtain  an  improved 
estimate  for  the  second  derivative .  The  process  can  be  repeated  until  a 
convergence  criterion  is  satisfied.  Related  expressions  can  be 
obtained  for  the  higher-order  derivatives  of  Eq.  (3);  however,  the 
expressions  are  quite  lengthy. 

The  computational  cost  of  evaluating  thermal  properties  at  local 
conditions  is  negligible,  requiring  two  calls  to  a  variable  properties 
function  per  time  step.  The  cost  of  including  temperature  gradients  is 
larger.  The  calculation  of  the  spatial  derivatives  of  temperature  points 
requires  the  evaluation  of  all  of  the  c,(f)  coefficients,  then  new 
coefficients  for  the  digital  filter  must  be  calculated  and,  finally,  the 
filter  calculation  must  be  repeated.  The  additional  steps  increase  the 
processing  time  by  approximately  a  factor  of  4.  In  most  cases,  the 
computational  time  will  still  be  negligible.  However,  it  is  possible  to 
include  a  flag  to  turn  on  the  gradient  calculations  only  when  heat  flux 
exceeds  the  critical  value. 

Example  calculations  were  performed  to  assess  the  effectiveness 
of  the  variable  properties  approximations.  The  reference  case  was  a 
square-wave  pulse  of  heat  of  duration  Ax\/a.  This  case  is  useful  for 
conveying  the  time  response  and  the  accuracy  of  the  methods,  and  it 
is  an  idealized  version  of  the  heat  flux  pattern  encountered  in  tests  of 
heat-sink  rocket  chambers.  Test  case  data  consisting  of  temperatures 
at  the  measurement  locations  were  generated  using  a  numerical 
solution  of  the  nonlinear  heat  equation  with  functions  for  the 
temperature-dependent  properties  of  copper.  Figure  1_4  includes 
results  for  the  P 5  model  for  three  methods  of  handling  thermal 
properties  for  a  heat  flux  of  107  W/m2,  and  Fig.  _15  has  the  same 
information  for  108  W/m2.  At  107  W/m2,  the  effect  of  property 
variations  for  the  short  duration  pulse  is  negligible,  and  all 
approximations  give  the  same  result.  The  rise  time  is  approximately 
0.03  s,  and  the  heat  flux  is  reproduced  accurately.  For  this  case,  the 
temperature  rise  at  the  surface  was  1 13  Kin  0.13  s.  At  108  W/m2,  the 
effect  of  variable  properties  is  observable.  The  constant  properties 
approximation  does  not  converge  to  the  correct  value  but  rises 
steadily  over  the  duration  of  the  pulse,  and  the  fall  time  is  sig¬ 
nificantly  increased.  The  local  properties  approximation  more  accu¬ 
rately  captures  the  steady  plateau,  and  the  fall  rate  is  improved. 
However,  the  heat  flux  is  underpredicted  by  2-5%.  The  gradient 
approximation  exhibits  the  same  rise  and  fall  rates  as  the  local 
properties  approximation,  and  the  accuracy  is  somewhat  improved. 
However,  the  square  profile  is  not  as  well  defined  as  it  was  for  a  heat 
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Fig.  14  Effect  of  temperature-dependent  properties  on  heat  flux 
calculation  at  107  W/m2. 
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Fig.  15  Effect  of  temperature-dependent  properties  on  heat  flux 
calculation  at  108  W/m2. 


Fig.  16  Comparison  with  function  specification  with  three  future  time 
steps  for  heat  flux  of  108  W/m2. 


flux  of  1 07  W/m2 .  For  this  case,  the  surface  temperature  increases  by 
approximately  1200  K  over  0.13  s. 

F.  Comparison  with  Beck’s  Function  Specification  Method 

In  this  section,  we  compare  the  polynomial  model  with  a  widely 
used  inverse  heat  transfer  method,  the  function  specification  with 
future  time  steps  method,  for  a  case  in  which  temperature-dependent 
properties  are  significant.  The  method  is  described  fully  in  [8].  The 
method  was  implemented  using  The  MathWorks  MATLAB© 
version  2008  library  functions.  The  pdepe  nonlinear  equation  solver 
was  used  to  model  the  heat  equation.  The  pdepe  function  solves 
initial-boundary  value  problems  for  systems  of  parabolic  and  elliptic 
partial  differential  equations  in  one  space  variable  [25],  The  spatial 
domain  was  divided  into  forward  and  inverse  regions,  and  the  inverse 
region  between  the  sensor  and  the  surface  consisted  of  five  finite 
volumes.  The  fzero  function  was  used  to  find  the  minimum  of  the 
convergence  criteria.  The  fzero  function  uses  a  combination  of 
bisection,  secant,  and  inverse  quadratic  interpolation  methods  [26]. 
The  size  and  number  of  future  time  steps  were  based  on  optimal 
values  contained  in  Sec.  5.6.1  of  [8],  For  the  cases  described  next,  the 
number  of  time  steps  was  four,  and  the  length  of  the  time  step  was 
QA5x\/a.  Figure  H)  contains  the  results  of  an  example  calculation 
for  a  heat  flux  of  108  W/m2.  The  function  specification  method 
accurately  reproduces  the  steady  portion  of  the  heat  flux  and  has 
similar  rise  and  fall  times  as  the  Ps  model.  The  function  specification 
method  achieves  greater  accuracy  by  modeling  the  full  domain.  The 
computational  cost  of  the  higher  accuracy  in  terms  of  run  time  was 
approximately  a  factor  of  105  on  this  example  problem. 

III.  Conclusions 

A  method  has  been  described  for  using  approximating 
polynomials  to  solve  the  inverse  heat  transfer  problem  for  the  case 
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of  two  temperature  sensors  embedded  in  the  wall  of  a  chamber.  The 
method  does  not  require  surface  junction  thermocouples,  which  are 
prone  to  failure  and  produce  noisy  signals  in  rocket  engine  flows,  and 
is  well  suited  for  studies  of  the  effects  of  surface  features  on  heat 
transfer  enhancement.  An  approximating  polynomial  is  constrained 
to  match  the  temperatures  and  the  even-numbered  spatial  derivatives 
at  the  measurement  points.  The  method  requires  only  current  values 
of  temperature  and  its  rate  of  change,  and  the  boundary  and  initial 
conditions  are  arbitrary.  The  algorithm  can  be  represented  as  a  low- 
pass  filter,  and  the  gain  and  phase  behavior  have  been  characterized. 
The  placement  of  the  sensors  affects  the  frequency  cutoff  and  the 
noise  response,  and  optimum  values  for  the  relative  positions  of  the 
sensors  have  been  obtained.  The  method  is  computationally  efficient, 
requiring  approximately  100  multiply  and  add  operations  per 
measurement,  enabling  real-time  processing  in  high-channel-count 
systems. 
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